\chapter{\label{chp:pmc}The Time-gated Perturbation Monte Carlo Method}
In this chapter we present the development and the first use of the time-gated perturbation Monte Carlo (pMC) method for 3D reconstruction of functional parameters in small animals. This method is parallelized to efficiently calculate the Jacobian utilizing the photon path information. The performance of this method is assessed \emph{in silico} by reconstruction of the distribution of intrinsic optical properties, and then by direct functional imaging with geometric and spectrum priors under different probe settings. Lastly, a time-resolved pattern illumination scheme is presented and the proposed method is experimentally validated\let\thefootnote\relax\footnotetext{
Portions of this chapter previously appeared as:
\begin{itemize}
 \item[] J. Chen and Xavier Intes, ``Time-gated perturbation Monte Carlo for whole body functional imaging in small animals,'' \textit{Opt. Express}, vol. 17, no. 22, pp. 19566-19579, 2009.
 \item[] J. Chen \textit{et al.}, ``Time-resolved diffuse optical tomography with patterned-light illumination and detection,'' \textit{Opt. Lett.}, vol. 35, no. 13, pp. 2121-2123, 2010.
 \item[] V. Venugopal, J. Chen and X. Intes,``Development of an optical imaging platform for functional imaging of small animals using wide-field excitation,''\textit{Biomed. Opt. Express}, vol. 1, no. 1, pp. 143-156, 2010.
\end{itemize}
}.

\section{\label{sec:meth_tr_mc} Introduction to the Monte Carlo method for light propagation}
% historical review
The Monte Carlo method originated at the Los Alamos Scientific Laboratory by Stanislaw Ulam, who suggested the simulation of neurons' flight paths by random sampling. John von Neumann, who later developed a detailed proposal for this method, named it Monte Carlo in reference of the Monte Carlo casino in Monaco \cite{Anderson1986}. Monte Carlo methods were first used for military purposes, such as the hydrogen bomb development in the Manhattan project during the 1950s. They then began to be applied in a wide range of applications in other fields of study, such as mathematics, physics, physical chemistry, and operations research. Monte Carlo methods had been employed for photon propagation in turbid media by pioneers from the late 1970s \cite{Witt1977,Prahl1989, Wilson1983}. Currently, a few Monte Carlo computation packages, including MCML by Lihong Wang et al. for multi-layered geometry \cite{Wang1995}, and tMCimg for 3D voxel-based time-gated light transport \cite{Boas2002}, are publicly available. Details of the light propagation model are well established and can be found in \cite{Prahl1989}. In this section, a brief review of basic concepts of probability theory applied in the Monte Carlo method is provided.

% probability theory and general MC
All Monte Carlo simulations involve generating random variables that follow certain probability distribution rules. To achieve this, inverse probability transform is usually used to convert the uniformly distributed random numbers, that can be easily obtained from random number generators, to the given probability distribution. This method is based on the fact that the cumulative distribution function $F_X$ of a continuous random variable $X$ has a uniform distribution on [0, 1]. If $F_X^{-1}$ is the inverse of $F_X$, a random variable $X' = F_X^{-1}(Y)$ with $Y$ having a uniform distribution on [0, 1], will have the same cumulative distribution function of the random variable $X$.

In Monte Carlo methods for photon propagation, the most important probability distribution resides in those describing the step size of photon movement between two scattering events, and the angle of deflection in a photon's trajectory after a scattering collision. Take the calculation of the step size as an example: from (\ref{eq:prob_scat_length}) the cumulative distribution of $P_{l_s}$ (probability density function of scattering step size) is given by:
\begin{equation}
F_{l_s}(l_s) = \int_0^{l_s}P_{l_s}(l_s)dl_s = 1-\exp(-\mu_s {l_s}).
\end{equation}
Replace $F_{l_s}(l_s)$ by a random variable $\xi=F_{l_s}(l_s)$ with uniform distribution between 0 and 1, we have
\begin{equation}
\label{eq:scat_l_xi}
\Delta s = -\frac{\ln{\xi}}{\mu_s},
\end{equation}
where $\Delta s$ is the calculated step size, a random variable that has a distribution identical with the actual scattering length $l_s$. Similarly, by using the Henyey-Greenstein phase function (\ref{eq:henyey}), the angle after a scattering event can be calculated as:
\begin{equation}
\label{eq:henyey_xi}
\cos\theta =\left\{
\begin{array}{ll}
 \frac{1}{2g}[1+g^2-(\frac{1-g^2}{1-g+2g\xi})] & \text{if } g \neq 0,\\
 1-2\xi &\text{if } g=0,
\end{array}
\right.
\end{equation}
where $\xi$ is again a random variable uniformly distributed between 0 and 1.

Fig. \ref{fig:MC_flow} summarizes the key steps in Monte Carlo simulations \cite{Prahl1989, Wang1995, Boas2002}. The initial position of a photon is set as the injection point of a punctual light beam and the initial time is set at $t_0=0$, to simulate the spatial and temporal impulse of a source. The initial weight of the photon is assigned as 1. As it moves, the distance to the next scattering event and the scattering angle are calculated by (\ref{eq:scat_l_xi}) and (\ref{eq:henyey_xi}), respectively. Upon reaching the internal or surface boundary, the photon is either reflected or refracted off of the interface. The probability and the angle of reflection or refraction are respectively given by the Fresnel equations and the Snell's law. As the photon passes from voxel to voxel, the photon weight is attenuated exponentially according to the local absorption coefficient and the distance traveled. The photon position and weight are recorded at regular time intervals to update the time-domain fluence. The photon is repeatedly moved until it either escapes from the boundary or its propagation time exceeds the maximum time of interest. Photons are continuously launched until the desired number of photons has been propagated.
\begin{figure}[htbp]
\centering
\includegraphics{pmc-6}
\caption{Flow chart of the time-gated Monte Carlo method. $\Delta s$ is the distance to the next scattering event.}
\label{fig:MC_flow}
\end{figure}

\section{\label{sec:meth_tr_mc} Using time gates in the perturbation Monte Carlo method}
The determination of spatially distributed optical parameters is important in near-infrared small animal optical imaging. Although investigators have determined the bulk optical properties based on homogeneous models \cite{Pan2007,Hayakawa2001}, direct imaging of the optical properties of complex heterogeneous tissue remains a difficult task.

% pMC history
The perturbation Monte Carlo (pMC) method is an analytical formulation for retrieving the optical properties based on Monte Carlo simulations. This method assumes the heterogeneous change over the background is sufficiently small that the derivative information of a background forward simulation can be used to estimate the perturbation. Early descriptions of obtaining the derivatives of forward Monte Carlo methods can be found for neutron transport applications, in works by Dejonghe \cite{Dejonghe1982} and Rief \cite{Rief1984} for sensitivity analysis, and Hall \cite{Hall1982} for adjusting material properties to fit reaction rates. For optics, pMC was first developed for bulk optical properties estimation in layered media by Hayakawa~\cite{Hayakawa2001} (cervical stromal like tissue), and was soon extended to 2D grid tomographic settings using simulated measurements for low-scattering media where the DE fails~\cite{Kumar2004}. This method has undergone experimental validation under these settings~\cite{Seo2007,Yalavarthy2005}, and has been applied to visualize brain hemodynamics related absorption changes recently~\cite{Heiskala2005}.

% pMC principles
The principle of pMC is summarized as follow: with the initial weight $w$ of one photon, the number of scattering events $P$ and the path length $L$ in a perturbed region (changes in optical properties), a prediction of the measurement of one photon at the detector can be estimated by:
\begin{equation}\label{eq:pMC_analytic}
{\hat w} = w\left( \frac{{\hat \mu } _s /{\hat \mu } _t }{{\mu _s /\mu _t }}
\right)^P \left( \frac{{\hat \mu } _t }{\mu _t } \right)^P \exp \left( { \left(
-{\hat \mu } _t- \mu _t \right)L} \right),
\end{equation}
where ${\hat \mu } _s = \mu _s + \delta \mu _s$, ${\hat \mu } _a = \mu _a + \delta \mu _a$, $\mu _t = \mu _a + \mu _s$ and $\hat{\mu} _t = \hat{\mu} _a + \hat{\mu} _s$ , with $\mu _a$ and $\mu _s$ being the initial values of the absorption and scattering coefficient and $\delta \mu_a$ and $\delta \mu_s$ the corresponding perturbations. This equation provides an accurate way to obtain the exit weight for every single photon. Based on the analytical solution (\ref{eq:pMC_analytic}), a simulation recording $w$, $P$ and $L$ for each photon with a given set of absorption and scattering coefficients, can be employed for optical property estimation in the perturbed regions, which is the fundamental of previous pMC \cite{Hayakawa2001,Kumar2004}.

Notice here if multiple perturbations are encountered, the predicted measurement $\hat{w}$ of this photon should be scaled by a product of all the perturbed regions. Each voxel can be considered as a perturbed region potentially, thus the number of unknowns is determined by the number of voxels in the region of interest, with the final estimation of the detected weight $\hat w$ written as:
\begin{equation}\label{eq:pMC_t_analytic2}
{\hat w} = w\prod_{j=1}^n\left( \frac{{\hat \mu } _s(r_j) /{\hat \mu }
_t(r_j) }{{\mu _s(r_j) /\mu _t(r_j) }} \right)^{P(r_j)} \left( \frac{{\hat
\mu } _t(r_j) }{\mu _t(r_j) } \right)^{P(r_j)} \exp \left( { \left( -{\hat
\mu } _t(r_j)- \mu _t(r_j) \right)L(r_j)} \right),
\end{equation}
where $P(r_j)$ is the number of collisions and $L(r_j)$ is the path length at voxel $r_j$ $(j = 1,...,n)$ for this photon.

It is straightforward to extend the pMC method to time-resolved studies. Assuming that the perturbations are time invariant over the propagation time (\~ns), the time dependent photon weight becomes:
\begin{equation}\label{eq:pMC_t_analytic1}
{\hat w(t)} = w(t)\prod_{j=1}^n\left( \frac{{\hat \mu } _s(r_j) /{\hat \mu }
_t(r_j) }{{\mu _s(r_j) /\mu _t(r_j) }} \right)^{P(r_j)} \left( \frac{{\hat
\mu } _t(r_j) }{\mu _t(r_j) } \right)^{P(r_j)} \exp \left( { \left( -{\hat
\mu } _t(r_j)- \mu _t(r_j) \right)L(r_j)} \right),
\end{equation}
where $\hat w(t)$ and $w(t)$ is the photon weight detected in a time-bin $t$.

The measurement at $t$ can thus be determined by summing up all $N$ photons detected at this time:
\begin{equation}\label{eq:pMC_t_analytic2}
\begin{split}
{\hat W(t)} = \sum_{K=1}^N w_K\prod_{j=1}^n\left( \frac{{\hat \mu } _s(r_j) /{\hat \mu }
_t(r_j) }{{\mu _s(r_j) /\mu _t(r_j) }} \right)^{P_K(r_j)} \left( \frac{{\hat
\mu } _t(r_j) }{\mu _t(r_j) } \right)^{P_K(r_j)}\\
\exp \left( { \left( -{\hat
\mu } _t(r_j)- \mu _t(r_j) \right)L_K(r_j)} \right),
\end{split}
\end{equation}
where the subscript $K = 1,...,N$ denotes the index of the photons and $\hat W(t)$ the measurement at $t$.

With $\delta \mu(r_j)$ the perturbation in optical coefficients (either $\mu_a$ or $\mu_s$) at $r_j$ , (\ref{eq:pMC_t_analytic2}) provides a means to estimate the changes in photon weight introduced by changes in $\delta \mu(r_j)$ and is then used to form the time-gated Jacobian matrix $J(t)$ with $J_{ij}(t)=\partial \hat W_i(t)/\partial\delta\mu(r_j)$ for this time gate. At a specific time-gate $t$ , the problem can be then casted in the formulation:
\begin{equation}\label{eq:linear1}
\left[\begin{array}{c}
 \delta W_1(sd,t) \\
 \vdots \\
 \delta W_m(sd,t)
 \end{array}\right]
 =
 \left[\begin{array}{ccc}
 J_{11} &\dots&J_{1n}\\
 \vdots &\ddots&\vdots\\
 J_{1m}&\dots&J_{mn}
 \end{array}\right]
 \times
 \left[\begin{array}{c}
 \delta\mu(r_1) \\
 \vdots \\
 \delta\mu(r_n)
 \end{array}\right]=J(t) \times
 \left[\begin{array}{c}
 \delta\mu(r_1) \\
 \vdots \\
 \delta\mu(r_n)
 \end{array}\right],
\end{equation}
where $\delta W_i(sd,t) = \hat W_i(sd,t) - W_i(sd,t)(i = 1,...,m)$, is the difference of unperturbed and perturbed measurement for source-detector pair $sd$ at $t$. This equation can be easily expanded for multiple gates:
\begin{equation}\label{eq:linear1}
\left[\begin{array}{c}
 \delta W_1(sd,t_1) \\
 \vdots \\
 \delta W_{m\times k}(sd,t_k)
 \end{array}\right]
 =
 \left[\begin{array}{c}
 J_{t_1}\\
 \vdots\\
 J_{t_k}
 \end{array}\right]
 \times
 \left[\begin{array}{c}
 \delta\mu(r_1) \\
 \vdots \\
 \delta\mu(r_n)
 \end{array}\right]=J \times
 \left[\begin{array}{c}
 \delta\mu(r_1) \\
 \vdots \\
 \delta\mu(r_n)
 \end{array}\right],
\end{equation}
where $k$ is the number of selected gates and $m\times k$ is the total number of selections of measurements.

For a computational experiment, two sets of Monte Carlo simulations are required to retrieve optical property maps according to the method discussed above. One is used to perform a simulation with an initial guess of (or background) optical properties, which does not have to be homogeneous; any heterogeneity can be considered. This simulation generates the unperturbed detector reading $W$ and the Jacobian $J$. The second set is to simulate the perturbed measurement ($\widehat{W} $) with changes in optical properties. This set of simulations has the benefit of providing measurements based on the most rigorous forward model of light propagation. In real experiments, the unperturbed and perturbed detector readings are replaced by the actual measurements and only the simulation for calculating Jacobian is required. A linear system can be formed after these two simulations are finished.

Some of the applications of pMC are based on a nonlinear optimization framework \cite{Kumar2004, Seo2007}, that is, to perform Monte carlo simulations using the optical property estimation at each iteration to get the derivatives. This method is not limited by the constraint that changes of heterogeneity have to be sufficiently small in the perturbation method. However, in our 3D, time-gated applications, Monte Carlo simulations for each iteration would take an extremely long time. To alleviate this drawback, the White Monte Carlo method can be applied herein for scaling the reconstructed absorption coefficients \cite{Alerstam2008a}. For each iteration step, the photon path information is first utilized to estimate perturbed maps of optical properties. Then, the perturbed optical properties and the saved photon path information are employed to calculate the perturbed weights by (\ref{eq:pMC_analytic}). These perturbed optical properties and perturbed weights are set to the initial optical properties and weights for the next iteration step. The iterations will stop when the difference between the estimated weights and the actual measurements is sufficiently small. However, it is not applicable for dramatic changes in scattering coefficients due to the changed photon paths.

To summarize, the above method is depicted in Fig.
\ref{fig:ab_s_flow}: based on the forward Monte Carlo simulation, the Jacobian is obtained; data processing techniques such as pre-conditioning of the matrix is used to reduce the ill-posedness of the inverse problem; finally the 3D distribution of optical properties can be retrieved by solving the linear system.
\begin{figure}[htbp]
\centering
\includegraphics{ab_s_flow}
\caption{Flowchart for the absorption and scattering reconstruction.}
\label{fig:ab_s_flow}
\end{figure}

% the need for parallelization
However, applying the time-gated pMC method to 3D tomography faces several challenges in its technical implementation. First, the size of the Jacobian in time-gated tomography with high-resolution 3D voxelized geometry is extremely large compared to that in basic pMC applications. This difference is mainly due to the dimension increase of the Jacobian in space and time compared to layered medium (1D) or grided maps (2D) in CW mode. The size of the matrix storing the photon paths and number of scattering events for reconstruction is dependent on the product of the voxels of the imaging volume, the number of detectors and the number of gates. Consider a 4cm (width)$\times$ 2cm (height) $\times$ 10cm (length) small animal model with 1mm$^3$ voxels and 100 gates, the size of the the Jacobian is increased by $100\times100=10,000$ times that in a 2D pixelized CW problem. In addition, more sources and detectors have to be employed to resolved 3D images, which results in further increase of Jacobian size by $\sim100$ folds.

Second, as mentioned previously, much more photons need to be simulated due to the decreased statistical reliability in time-resolved Monte Carlo simulations, especially for generating high-resolution temporal profiles. This not only translates into larger computational time, but also precludes the direct storage of the photon path length and scattering events for each photon, which is essential in pMC method.

Due to the reasons outlined above, storing each photon path voxel-wise and sort them for different gates and source-detector pairs, as in the classic pMC applications (as appeared in (\ref{eq:pMC_t_analytic2})), becomes infeasible. To overcome this problem, the accumulated detector readings and the weighted average of the collected photon information is stored at each voxel for each time-gate and source-detector pair. Equation (\ref{eq:pMC_t_analytic2}) is thus modified to:
\begin{equation}\label{eq:pMC_t_analytic3}
\begin{split}
{\hat W(t)} = W(t)\prod_{j=1}^n\left( \frac{{\hat \mu } _s(r_j) /{\hat \mu }
_t(r_j) }{{\mu _s(r_j) /\mu _t(r_j) }} \right)^{\overline{P}(r_j,t)} \left( \frac{{\hat
\mu } _t(r_j) }{\mu _t(r_j) } \right)^{\overline{P}(r_j,t)}\\
 \exp \left( { \left( -{\hat
\mu } _t(r_j)- \mu _t(r_j) \right)\overline{L}(r_j,t)} \right),
\end{split}
\end{equation}
where $\overline{P}(r_j,t)$ and $\overline{L}(r_j,t)$ are the time-dependent weighted average of the scattering and path length information. However, this matrix needs to be kept in memory during the simulation, and the memory capacity on a single computer is not sufficient to store it. Based on this consideration, together with the goal to calculate pMC-based Jacobian in a computational efficient fashion, a parallelization scheme is employed.

\section{Parallelization and computational feasibility}
\subsection{Parallelization scheme}
This parallelization scheme is two-fold. First, the large number of samplings in Monte Carlo algorithms can be easily divided into smaller groups for separate calculation on different CPUs (computation nodes). This makes forward Monte Carlo simulations of light propagation ideally suited for parallel computing by assigning only a portion of the photons to each CPU, which has been widely used \cite{Boas2002}. To assure the randomness, each CPU has a distinct seed for the random number generator for different pseudo-random sequences. The number of photons that each node calculates is obtained by dividing the total number of photons by the total number of nodes, and assigned at the beginning of computation. When this number is sufficiently large ($>100$), which is the case in our application, the computational imbalance between nodes is negligible due to the randomness of photon paths. This method circumvents the difficult implementation of dynamic assignment while achieving similar parallelization efficiency.

Second, the matrix of photon path length and scattering information is allocated distributively, which is an unique feature that sets this pMC code apart from the other Monte Carlo codes. A few ``master nodes'' are assigned for the purpose of storing these matrices, with each master node containing time-gated matrices from one or more detectors, depending on the problem size. Upon allocation of master nodes, other ``slave nodes'' begin the simulation of photon propagation. If a photon arrives at a detector, the photon path and scattering information are stored in a buffer with a memory-efficient data structure, due to the sparsity of a 3D voxel-based matrix of a photon's path. Once the buffer is full, this information and the photon's arrival time is subsequently transferred to the corresponding master node based on the detector it reaches. Master nodes are responsible for converting and adding the data received to the matrix, and output all the matrices to the disk when calculation is complete.

This parallelization is implemented using the Message Passing Interface (MPI) protocol for C. MPI is an application programming interface that provides basic data transferring and flow control function for both massively parallel machines and on workstation clusters. Communication functions MPI\_send() and MPI\_recv() are called for passing the photon data to master nodes, and MPI\_barrier to make sure the same pace was present on all nodes. MPI\_reduce() is also called for adding up the results from different nodes. The parallel implementation scheme is provided in Fig. \ref{fig:paral_pmc}.
\begin{figure}[htbp]
\centering
\includegraphics{pmc-7}
\caption{Implementation of the parallel pMC code.}
\label{fig:paral_pmc}
\end{figure}

\subsection{Computational efficiency}
We evaluated the computational efficiency of the time-gated forward MC versus the number of CPUs, number of gates and number of detectors used in the simulations. The motivation of
this evaluation is due to (a) the low computational efficiency of MC method and (b) the temporally and spatially dense datasets required in time-resolved whole-body imaging.

First, we tested the computational feasibility of the proposed method using different numbers of CPUs for one source-detector pair. In this case, the source and detector were facing each other in a transmittance geometry (2cm separation), probing an homogenous media. Second, we evaluated the effect on time cost caused by changes in the number of detectors and gates simulated. The recent development in the ultra-fast gated CCD camera makes acquiring spatially dense time-resolved datasets possible. Such sub-millimeter spaced arrays of sources and detectors (that is, data to the order of 10$^{4}$ - 10$^{6}$ measurements or more) are required for high-fidelity, small animal imaging \cite{Konecky2008}. In the pMC simulations, multiple detector readings and the corresponding Jacobian can be obtained simultaneously for all gates. We investigated the time cost associated with simulating 16, 64, 256 and 1024 evenly spaced detectors. Sources and detectors were simulated with a 1mm diameter in all simulations herein. The area covered by different detector sets was fixed to 16mm $\times$16mm. Additionally, for each of these source-detector settings, time-gated Jacobians for 1, 10, 20, 30, 40, 50 gates were recorded and the time cost for each simulation was evaluated.

\begin{figure}[htbp]
 \centering
 \subfigure[ ]{
 \label{fig:time_costa}
\includegraphics[width=7cm]{time_cost}}
 \subfigure[ ]{
 \label{fig:time-costb}
 \includegraphics[width=7cm]{time_nodes}}
 \subfigure[ ]{
 \label{fig:time_det_gate}
 \includegraphics[width=7cm]{time_effi_det}}
 \caption{Time efficiency estimation. (a) Time cost using increasing number of photons (from $10^{5}$ to $10^{10}$) with homogeneous optical properties on a 64-CPU cluster (Borg at SCOREC in RPI); (b) Time cost with increasing number of nodes (from 512 to 4096); (c)The computational time simulating different gates and detectors.}
 \label{fig:time_efficency}
\end{figure}

Fig. \ref{fig:time_efficency}a shows the time cost of simulations using increasing number of photons. When the simulated number of photons increases, the calculation time increases correspondingly. Fig. \ref{fig:time_efficency}b depicts the time cost using different number of nodes. On a single CPU, a  time-resolved Monte Carlo simulation for one source and one detector takes 62 hours to complete while it can be completed in under 5 minutes if 1024 nodes are used on BlueGene. It clearly shows the significant increase of computational efficiency thanks to parallel computing techniques.

Fig. \ref{fig:time_efficency}c shows the time cost of simulations using different numbers of detectors and gates. It is worthwhile to note that the time cost increases non-linearly as the gate number increases, which indicates that the proposed approach allows one to simulate large number of gates without significant penalty. Moreover, although parallel data transfer occurs in parallel computing, the time shift of using different numbers of detectors in the forward model is within 10\%, which allows simulations of dense detectors simultaneously. These two phenomena highlight the potential of time-gated pMC modeling for processing dense spatial and temporal datasets such as those obtained with a time-gated CCD camera.

\section{Reconstructing intrinsic optical parameters}
\subsection{Reconstructing single parameters \label{sec:pmc-single}}
In this section we describe the reconstruction of a single optical property perturbation (either the absorption or scattering coefficient) on a synthetic mouse model. \\

\noindent {\bf Computational settings}

A 3D mouse atlas created from CT and cryosection slices at University of Southern California \cite{Dogdas2007,Stout2002} was used to generate the anatomically accurate model to perform the MC simulations. The mouse model was discretized to $91\times35\times22$ voxels with size 1mm$^3$. The optical properties were considered constant in each voxel and may be different from other voxels.

In the unperturbed simulation, the whole mouse was considered homogenous, and each voxel had the same values $\mu_a = 0.04$cm$^{-1}, \mu_s = 40$cm$^{-1}$. We limited the heterogeneity to the lungs to calculate the perturbed measurements. Only one type of inhomogeneity (either $\mu_a=0.4$cm$^{-1}$ or $\mu_s=60$cm$^{-1}$) in lungs was considered at a time. The anisotropy factor $g$ was kept constant at 0.90 over the full body and the refractive index of the tissue at 1.37. The reduced scattering coefficients $\mu_s'$ are relatively low in this simulation compared to the average value in small animals but allow for faster computation (analyzed in Chapter \ref{chp:comp}).

In this example, a transmittance geometry with conventional point sources and detectors was employed. An array of 25 impulse sources on the mouse surface and 25 detectors with a 1mm radius are shown in Fig. \ref{fig:pmc-mouse-model}. The source-to-source and detector-to-detector distances were 3mm.\\

\begin{figure}[htbp]
 \centering
 \includegraphics{pmc-8}
 \caption{(a) The mouse model with the extracted skin-shape and lungs used in this work; (b) source-detector geometry; (c) time-resolved measurements obtained from Monte Carlo simulations at the center source-detector pair.}
 \label{fig:pmc-mouse-model}
\end{figure}

\noindent {\bf Data processing}
Following the same noise type of real light propagation, the noise in Monte Carlo simulations is characterized as a Poison noise \cite{Dereniak1984,Boas2001}. The presence of noise degrades the pMC-based Jacobians' spatial integrity and makes the reconstruction less stable, due to the ill-posedness of the inverse problem. To obtain stability in reconstruction algorithms, the selection of source-detector pairs and time gates with meaningful statistics is required. The most effective method to achieve this is using a photon count threshold according to the TPSF to rule out the gates with low SNR. Moreover, in this reconstruction, the Anscombe transformation was used to transform the Poison noise into an approximate Gaussian additive noise \cite{Anscombe1948,suzuki1985}, based on:
\begin{equation}
b = 2\sqrt {a + \frac{3}{8}},
\end{equation}
where $a$ is the variable with Poisson distribution and $b$ is a Gaussian additive white noise, which can be reduced by many statistical processing tools. Previous experimental results indicates the use of Wiener filter with Anscombe transformation gives superior results \cite{Mascarenhas1999}, therefore we used the 3D Wiener Filtering procedure to pre-process the Jacobian. The inverse Anscombe transformation was then performed to obtain the final Jacobians to be used in the inverse problem.\\

\noindent {\bf Solving the inverse problem}

The derivative information, denoted as the Jacobian $J$, of the signal summation $W$ at any detector can be drawn from the analytical expression in (\ref{eq:pMC_t_analytic3}). When $J$ is not square (typically $n<<m$ in DOT), a method seeking for a solution can be employed by minimizing an objective function $\Psi$, commonly in a least square form:
\begin{equation}
\label{fig:cg} \Psi = ||\Delta W-J\Delta \mu||,
\end{equation}
where $\Delta W$ is the difference between the measured signal $\hat{W}$ and the data predicted by the forward model using (\ref{eq:pMC_analytic}) on the initial estimate $(\mu_a,\mu_s)$. To solve this optimization problem, the conjugate gradients method for normal equations, which is ideal for large-scale linear systems due to its iterative nature, was employed. The algorithm stops if either 50 iterations was reached or the relative residue was below 0.01.\\

\noindent {\bf Statistical stability and denoising technique}

To estimate the required photons for statistically reliable results, simulations using $10^5-10^9$ photons for Jacobian calculation were casted into the reconstruction scheme to obtain the optical perturbation. All these reconstructions employed the CW measurements from simulations using $10^9$ photons to avoid the bias caused from measurement differences. Fig. \ref{fig:pMC-recon-stat} shows one slice of the reconstructed absorption or scattering coefficients, in which the perturbation only occurred in $\mu_a$ or $\mu_s$ using different number of photons. As expected, reconstructions with less noise can be observed when the number of simulated photons increases with minor improvements after $10^9$ photons, while the reconstruction with $10^8$ photons can already provide meaningful results.
\begin{figure}[htbp]
\centering
\includegraphics[width=\textwidth]{CW_slices}
\caption{Reconstructed absorption (top row) coefficients and scattering coefficients (bottom row) at slice $z=13$.}
\label{fig:pMC-recon-stat}
\end{figure}

The noise is thus a critical factor for the quality of reconstructions when less number of photons are launched. The impact of denoising technique on the Jacobian and reconstruction is shown in Fig. \ref{fig:anscombe}, where the maximum reconstructed values for absorption and scattering are improved by 19.3\% and 63.4\%, respectively. This implies similar reconstruction results can be achieved with filtered data using less number of photons thus shorter computational time.
\begin{figure}[htbp]
\centering\
 \includegraphics{pmc-9}
 \caption{ (a) Unfiltered and (b) filtered Jacobians for one source-detector pair and (c) the reconstructions for absorption (top row) and scattering (bottom row) using these Jacobians along the line $x = 11,z = 16$ ($10^{9}$ photons used).}
 \label{fig:anscombe}
\end{figure}\\

\noindent {\bf Time-gated reconstruction}

Based on the above results, we investigated the potential of reconstructing the optical properties using the time-gated data. The filtered time-gated Jacobians with $10^9$ photons at the rising half maximum, maximum and decaying half maximum of the TPSF were employed. The corresponding reconstructed images employing one time gate are presented in Fig. \ref{fig:tg_results}a, \ref{fig:tg_results}b and \ref{fig:tg_results}c, with the maximum reconstructed values 95.5\%, 103.8\%, 93.6\% for absorption and 97.4\%, 99.4\%, 92.3\% for scattering. Combining these three Jacobians provides qualitatively accurate result (Fig. \ref{fig:tg_results}d, quantification 97.9\% for absorption and 99.2\% for scattering) while retaining the resolution similar to the early-gate based reconstruction. The reconstructed iso-volumes are shown in Fig. \ref{fig:iso}, in which a structural similarity to the expected organs is observed.
\begin{figure}[htbp]
\centering\includegraphics[width=\textwidth]{pmc-10}
\caption{Reconstructed absorption (top row) and scattering (bottom row) images employing the Jacobians at (a) rising 50\%, (b) maximum, (c) decaying 50\% gates of the TPSF and (d) combined Jacobian of (a), (b) and (c).}
\label{fig:tg_results}
\end{figure}
\begin{figure}[htbp]
\centering\ \subfigure[ ]{
 \label{fig:iso_ab}
 \includegraphics[width=7cm]{recon_iso_ab}}
 \subfigure[ ]{
 \label{fig:iso_sc}
 \includegraphics[width=7cm]{recon_iso_sc}}
 \caption{Reconstructed (a) absorption and (b) scattering coefficients using three time gates. The edge of the inclusions is 1/3 of the reconstructed maximum value.}
 \label{fig:iso}
\end{figure}

As seen in Fig. \ref{fig:tg_results}, the time-gating techniques can provide additional information -- early photons are associated with structural information (scattering) and later photons with functional information (absorption). This information makes casting optical inverse problem within the time-gated framework highly beneficial. The combined gates yield more stable reconstructions to provide both anatomical and physiological information, and overcome the limitations in typical DE-based DOT systems which suffer from low resolution and partial volume effect.

\subsection{Simultaneous reconstruction of both absorption and scattering coefficients}
In preclinical scenarios, both absorbing and scattering heterogeneities may exist simultaneously, thus it is desirable to reconstruct the two parameters in the same inverse problem. The experimental acquisition using continuous wave systems is limited in providing absolute estimates and/or differentiate the tissue absorption and scattering properties \cite{Arridge1998}. Since time-domain technology offers a rich dataset that is unmatched by CW and frequency domain methods, it is the method of choice for reconstructing both parameters simultaneously.

In this section we investigated the performance of time-gated DOT in the presence of absorbing and scattering heterogeneities. DOT is by nature an ill-posed, ill-conditioned inverse problem, thus suffers from low resolution, partial volume effect and inter-parameter cross-talk. Combining measurements at several time gates in the inverse problem will allow for reconstructing both parameters simultaneously, where its ability to separate the optical parameters may be associated with the selection of appropriate time gates.

When we have variations in both absorption and scattering, the final weight of photons are perturbed by both changes. This leads us to the following matrix to invert:
\begin{equation}
\left[\begin{array}{c}
 \delta W_1(sd,t) \\
 \vdots \\
 \delta W_m(sd,t)
 \end{array}\right]
 =
 \left[\begin{array}{cccccc}
 J_{a11} &\dots&J_{a1n}& J_{s11} &\dots&J_{s1n}\\
 \vdots &&\ddots&&\vdots\\
 J_{s1m}&\dots&J_{smn} & J_{s1m}&\dots&J_{smn}
 \end{array}\right]
 \times
 \left[\begin{array}{c}
 \delta\mu_a(r_1) \\
 \vdots \\
 \delta\mu_a(r_n)\\
 \delta\mu_s(r_1) \\
 \vdots \\
 \delta\mu_s(r_n)
 \end{array}\right].
\end{equation}

%\subsubsection{\label{sec:sim_absc}simulation results}
The mouse model used in this investigation was the same as described in Section \ref{sec:pmc-single}. The background optical values are $\mu_a$ = 0.04cm$^{-1}$, $\mu_s'$ = 4mm$^{-1}$ with the anisotropy factor $g$ = 0.90 and the refractive index 1.37. The illumination/detector setting is depicted in Fig. \ref{fig:ab_s_setting}. Two heterogeneous inclusions were considered. One was absorptive (red) with a contrast of 10 ($\delta \mu_a$= 0.36cm$^{-1}$) and the other was scattering (blue) with a contrast of 2 ($\delta \mu_s'$ = 4cm$^{-1}$). We simulated an array of 25 point impulse sources and 25 detectors on the mouse surface with a 1mm radius and 2.5mm separation evenly spanned in the perturbed area, with 10$^{10}$ photons launched per source. Simulations were performed in transmittance geometry with 200ps gate-width and 20ps time steps. Three gates that corresponded to the half-maximum rising, maximum and half-maximum decaying gates as shown in Fig. \ref{fig:ab_s_recon}a were employed in the reconstructions. A 29mm $\times$ 29mm $\times$ 22mm area was reconstructed using a Conjugate Gradient algorithm with spatially variant regularization \cite{Pogue1999}. The reconstruction program was stopped at 50 iterations or when the relative error was below 0.01.
\begin{figure}[htbp]
\centering
\includegraphics[width=2.5in]{ab_s_setting}
\caption{Perturbative inclusions (blue: scattering; red: absorption), with sources and detectors geometry.}
\label{fig:ab_s_setting}
\end{figure}

As a typical ill-posed problem, the condition of the inverse problem is worsened when reconstructing both optical parameters. Preconditioning techniques can be applied in this case to balance the sensitivity of absorption and scattering by rescaling the weight matrix to be more uniform thus better-conditioned \cite{Pei2001}. In this study we employed column scaling, to normalize each column of the weight matrix to its average value \cite{Intes2005a}. This scheme has been proven to be the most effective scheme for providing a well-conditioned system.

To qualitatively assess the reconstruction quality, we defined the quantification of reconstructed $\mu_a$ as:
\begin{equation}
Q_{\mu_a}=\frac{\max[\Delta \mu_a(\Omega_{\Delta \mu_a})]}{H_{\Delta \mu_a}(\Omega_{\Delta \mu_a})}
 \times 100\%,
\end{equation}
and the crosstalk for ${\mu_a}$ as:
\begin{equation}
X_{\mu_a} = \frac{\max[{\Delta \mu_s}(\Omega_{\Delta \mu_a})]/H_{\Delta \mu_s}(\Omega_{\Delta \mu_s})}
 {Q_{\mu_a}}
 \times 100\%,
\end{equation}
where $\Omega_{\mu_a}$ denotes the known location of ${\mu_a}$ perturbation and $H_{\mu_a}$ denotes the expected value of ${\mu_a}$. The quantification and crosstalk for the $\mu_s$ was similarly evaluated.

Optical reconstructions based on time-resolved datasets are provided in Fig. \ref{fig:ab_s_recon} as well as CW reconstruction for comparison. The quantification and crosstalk using time-gated and CW data are summarized in Table \ref{tab:pmc-both1}. The enhanced performance of TG datasets over CW dataset is demonstrated in terms of quantification. For example, there is a 26\% error of the maximum reconstructed scattering value using time-gated data, compared to a 56\% error in the CW reconstruction. Moreover, the crosstalk comparison proves that the different types of heterogeneities were separated using the time-gated approach while the CW reconstruction was not able to identify between absorption and scattering.
\begin{figure}[htbp]
\centering
\includegraphics{ab_s_recon}
\caption{(a) Time response at the center detector and the selected gates. (b) Reconstructed $\Delta \mu$ map using time-domain data and CW data (showing 50\% to maximum only). Left column: time-domain reconstruction; right column: CW reconstruction; upper row: absorption; lower row: scattering.}
\label{fig:ab_s_recon}
\end{figure}

\begin{table}[htbp]
\caption{\label{tab:pmc-both1}{Quantitative comparison for simultaneous $\mu_a$ and $\mu_s$ reconstructions using TG and CW data.}}
\centering
\def\temptablewidth{0.7\textwidth}
{\rule{\temptablewidth}{1pt}}
\begin{tabular*}{\temptablewidth}{@{\extracolsep{\fill}}ccccc}
\multirow{2}{*}{Datatype} & \multicolumn{2}{c}{Quantif. (\%)}&\multicolumn{2}{c}{Crosstalk (\%)} \\
&${\mu_a}$&${\mu_s}$&${\mu_a}$&${\mu_s}$\\
\hline
TG & 93.6& 73.2& 8.4 & 13.6 \\
CW & 69.8& 43.7& 33.9& 173.9
\end{tabular*}
{\rule{\temptablewidth}{1pt}}
\end{table}

\section{Functional imaging using \emph{a priori} information}
In the previous sections the use of time-gated Monte Carlo method to reconstruct the optical coefficients is demonstrated. It is natural to expand this method to directly retrieve the functional parameters, such as the hemoglobin concentration. In this section, we consider a combination of the time-gated pMC method and an anatomy guided inverse problem formulation to achieve high resolution whole-body quantitative functional imaging for small animals. The direct multispectral functional imaging method and a hierarchical Bayesian reconstruction are applied on an anatomically accurate mouse model, as well as an empirical probe selection approach to improve efficiency and accuracy of reconstruction. The tomographic results using simulated data for the proposed approach is reported.

%%-----------------------------------------------------------
\subsection{Applying spectral \textit{prior} for direct functional imaging}

Assuming that deoxy-hemoglobin(Hb), oxy-hemoglobin(HbO$_2$) and water are the only chromophores contributing to the global absorption in the spectral window considered, the physiological information or the changes in the concentration of Hb, HbO$_2$ and H$_2$O at different wavelength for $j^{th}$ voxel are linearly related to the absorption coefficients as follow:
\begin{equation}\label{linear2}
 \delta\mu_a^{\lambda_k}(r_j)
 =
 \left[\begin{array}{ccc}
 \varepsilon_{Hb}^{\lambda_k} &\varepsilon_{HbO_2}^{\lambda_k}&\varepsilon_{H_2O}^{\lambda_k}\\
 \end{array}\right]
 \times
 \left[
 \begin{array}{c}
 \delta[ Hb(r_j)] \\
 \delta [HbO_2(r_j)] \\
 \delta [H_2O(r_j)]
 \end{array}\right],
\end{equation}
where $\delta \mu^{\lambda_k} _a (r_j )$ is the differential absorption coefficient at the $k^{th}$ wavelength $\lambda_k$ and for the $j^{th}$ voxel, $\varepsilon_{Hb}^{\lambda_i}$, $\varepsilon_{HbO_2}^{\lambda_i}$ and $\varepsilon_{H_2O}^{\lambda_i} $ are the extinction coefficient of the corresponding absorbing substances at the $k^{th}$ wavelength. Combining this relationship with (\ref{eq:linear1}), we obtain:
\begin{equation}\label{linear3}
\left[\begin{array}{c}
 \delta W_1^{\lambda_k}(sd,t) \\
 \vdots \\
 \delta W_m^{\lambda_k}(sd,t)
 \end{array}\right]
 =
 \left[\begin{array}{ccc}
 \varepsilon_{Hb}^{\lambda_k}J^{\lambda_k} &\varepsilon_{HbO_2}^{\lambda1}J^{\lambda_k}&\varepsilon_{H_2O}^{\lambda1}J^{\lambda_k}
 \end{array}
 \right]\times
 \left[\begin{array}{c}
 \delta[ Hb] \\
 \delta [HbO_2] \\
 \delta [H_2O],
 \end{array}\right]\end{equation}
where $\delta[ Hb] = [\delta[ Hb(r_1)] \dots \delta[ Hb(r_n)]]'$, $\delta[ HbO_2] = [\delta[ HbO_2(r_1)] \dots \delta[ HbO_2(r_n)]]'$ and $\delta[ H_2O] = [\delta[ H_2O(r_1)] \dots \delta[ H_2O(r_n)]]'$ are the perturbations in the chromophore concentration. $\delta W_1^{\lambda_k}(sd,t)$ and $J^{\lambda_k}$ are the differential detector reading and the Jacobian at the $k^{th}$ wavelength $\lambda_k$, respectively. Multispectral priors can thus
be incorporated into the tomographic inverse problem to directly retrieve functional maps \cite{Corlu2005}.

These parameters provide a direct way to access the blood volume (BV) and relative saturation of oxygen (StO$_2$) in tissue, which are BV = [Hb] + [HbO$_2$], StO$_2$ = HbO$_2$/BV.

% Jacobian calculation description
%%-----------------------------------------------------------
\subsection{Bayesian method for applying geometry \emph{prior}}
% Inverse iterative method description
The Bayesian formulation is used to incorporate anatomical and functional priors for overcoming the ill-posedness associated with whole-body imaging. This formulation constrains the space of unknowns with a spatially varying probability density function derived from high resolution anatomical maps and physiological priors that are implicitly defined to constrain the reconstructions in physiologically meaningful ranges. The full derivation of the theoretical developments for the inverse formulation and algorithm is presented in detail in references \cite{Guven2005, Intes2004}. A brief description of the method is provided below.

In order to estimate the hyperparameters, we adapt the conjugate gradient (CG) algorithm to include a hyperparameter estimation step followed by an image update. In this context, we apply an iterative empirical Bayesian approach to estimate the hyperparameters, which in turn gives the maximum a posteriori (MAP) estimates of the hyperparameters at each CG iteration prior to the image update:
\begin{equation}
\hat{x}_{MAP} = \arg \max_x\{\log p(y|x)+\log p(x, \sigma|C)\},
\end{equation}
where $y$ is the measurement vector, $p(y|x)$ is the data likelihood, $\sigma$ is the vector of variance of sub-images and $p(x,\sigma|C)$ is the conditional probability prior on the image $x$ and the hyperparameter $\sigma$ given the tissue label information $C$. The label information $C$ can be taken into account by considering one sub-image with the same label at a time. Since $p(x,\sigma|C) = p(x|\sigma, C)p(\sigma|C)$, the optimization problem can be solved by alternatively updating these two parameters image $x$
and the variance $\sigma$. Assuming a Gaussian distribution for $x$, the probability density function of $i^{th}$ sub-image $x_i$ is:
\begin{equation}
p(x_i,\sigma_i) = \frac{1}{(2\pi\sigma_i^2)^{N_i/2}}
\exp(-\frac{1}{2\sigma_i^2}\|x_i-\mu_i\|^2), i = 1,2,\ldots,M,
\label{eq:baye_2}
\end{equation}
where $M$ is the number of sub-regions, $N_i$ is the number of voxels in the $i^{th}$ sub-image, $\mu$ is the mean value of this sub-image. Similarly, assuming a Gaussian distribution for the
standard deviation $\sigma_i$ of the $i^{th}$ sub-image, we have:
\begin{equation}
p(x_i,\sigma_i) = \frac{1}{(2\pi\gamma_i^2)^{N_i/2}}
\exp(-\frac{1}{2\gamma_i^2}\|\sigma_i-\bar \sigma_i\|^2), i = 1,2,\ldots,M,
\label{eq:baye_3}
\end{equation}
where $\gamma_i$ is the variance and $\bar \sigma_i$ is the mean value of variance $\sigma_i$. $\gamma_i$ and $\bar \sigma_i$ are predefined as the physiological prior to constrain the reconstruction process. The image and noise are accommodated by respectively applying (\ref{eq:baye_2}) and (\ref{eq:baye_3}) for all sub-images at each update
along with the solution process.

% what is this????????????aaaaaaaaaaaaaaaaaaaa
A comparison of the reconstruction schemes for the optical properties and the functional parameters are depicted in Fig. \ref{fig:func_flow}. The main difference between them is the incorporation of the spatial and spectrum \textit{a priori} information for reducing the non-uniqueness of the inverse problem.
\begin{figure}[htbp]
 \centering
 \includegraphics{func_flow}
 \caption{\label{fig:func_flow} The flow chart for direct functional imaging compared to the reconstruction of the optical properties.}
\end{figure}

\subsection{\emph{In silico} investigation of geometry driven reconstruction} \label{config}
\subsubsection{Geometric model and probe settings}
We tested this pMC scheme using the anatomically accurate 3D mouse atlas described in Section \ref{sec:pmc-single}. To assess the performances of time-gated functional imaging based on pMC formulation for whole-body small animal imaging, we generated a synthetic mouse model bearing five organs (heart, stomach, liver, kidneys and non-inflated lungs). The original high-resolution and derived synthetic phantom employed in this work are provided in Fig. \ref{fig:mouse_model}. We limited the volume imaged to the upper torso-mid section (including kidneys), leading to a 33$\times$41$\times$20 voxels discretized volume to reconstruct.
\begin{figure}[htbp]
 \centering
 \includegraphics[width=\textwidth]{func_mouse_setting}
 \caption{(a) The original high resolution mouse model (b) The resized mouse model with five organs: heart (teal), stomach (green), liver (yellow), kidneys (red) and lungs (purple) extracted from (a). \label{fig:mouse_model}}
\end{figure}

Simulations were performed under a transmittance geometry mimicking the configuration of our time-resolved multi-spectral non-contact imaging platform \cite{Venugopal2009}. Due to memory limitations associated with the significant size of the linear system to invert, we selected a subset of 100 sources and 144 detectors for all reconstructions herein. We investigated two heuristic strategies for the selection of sources and detectors positions. First, we selected a linear distribution of sources and detectors that covered all the volume to be imaged. In this configuration, sources were evenly spaced at 4mm intervals whereas the detectors were evenly separated by 3.6mm (Fig. \ref{recon_nonlinear}a). Such homogenous distribution is the optimal configuration according to Singular Value Analysis (SVA)~\cite{Intes2005a}.

Second, the same number of sources and detectors were selected but with a non-linear spatial sampling distribution. Guided by the structural prior, denser source-detector pairs were employed around the area of significant overlap between the (i) lungs and the heart and (ii) the liver and the stomach. The denser source-to-source separation was 2mm and detector-to-detector separation was 1.8mm while in other area the separation of sources or detectors was doubled (Fig. \ref{recon_nonlinear}b).

\subsubsection{Functional and optical parameters}
The \emph{in vivo} functional properties of murine organs are still unknown to date. Average optical properties of small animals have been reported in transmittance but \emph{in vivo} organ
specific values have not been reported yet, due to the difficulty of resolving such parameters locally. To model realistic optical properties, we simulated optical properties based on
functional parameters derived from the literature ~\cite{Cheong1990,Alexandrakis2005}. Note that these values were obtained mainly from \textit{ex vivo} measurements and we expect significant deviations from \emph{in vivo} values.

Although the initial guess of optical properties can be heterogeneous in the unperturbed simulations, we considered homogeneous initial optical properties in the unperturbed simulations. That is, all organs are assigned identical optical parameters as the background.

Simulations were performed at six wavelengths evenly spanned from 700nm to 950nm. For simplicity, the scattering parameter was set to be spatially homogenous in both unperturbed and perturbed simulations, but to spectrally follow the power law in (\ref{eq:power_law}) with A = 10,600 and b =1.43 \cite{Cerussi2001}. The functional parameters that were simulated for each organs are provided in Table \ref{func_para}. These functional parameters were used to calculate the corresponding optical absorption values of each organ at the six wavelengths simulated, as tabulated in Table \ref{tab:pmc_opti_para}. For the volume considered herein and with the transmittance geometry, no difference was found between simulations based on the anisotropy factor $g = 0.9$ or $g = 0$. Thus g was kept constant at 0 over the full body to reduce the computational burden in the Monte Carlo simulations. The refractive index of the tissue was set at 1.37.

\begin{table}[htbp]
 \tabcolsep 0pt \caption{Functional parameters from literature.} \vspace*{-12pt}
\begin{center}
\def\temptablewidth{0.8\textwidth}
{\rule{\temptablewidth}{1pt}}
\begin{tabular*}{\temptablewidth}{@{\extracolsep{\fill}}ccccccc}
 &Background & Heart& Stomach & Liver & Kidneys &Lungs\\ \hline
 BV[$\mu$M] & 0.06 & 0.15 & 0.1 & 0.20 & 0.056 & 0.1\\
 StO$_2$(\%)&75\% & 75\%& 50\%&75\%& 75\%&85\%\\
 H$_2$O(\%)& 50&50&80&70&80&85
 \end{tabular*}
 {\rule{\temptablewidth}{1pt}}
 \end{center}
 \label{func_para}
 \end{table}

\begin{table}[htbp]
 \tabcolsep 0pt \caption{Absorption coefficients at different wavelengths used in the functional tomographic reconstruction.} \vspace*{-12pt}
\begin{center}
\def\temptablewidth{0.8\textwidth}
{\rule{\temptablewidth}{1pt}}
\begin{tabular*}{\temptablewidth}{@{\extracolsep{\fill}}ccccccc}
 $\mu_a$(cm$^{-1}$) &Background & Heart& Stomach & Liver & Kidneys &Lungs\\
\hline
 700nm & 0.217 & 0.538& 0.563& 0.718& 0.205& 0.281\\
 750nm&0.251 &0.607 & 0.536& 0.811& 0.243& 0.371\\
 800nm& 0.268& 0.655& 0.438& 0.874& 0.257& 0.450\\
 850nm & 0.332& 0.798& 0.501& 1.065& 0.324& 0.574\\
 900nm & 0.384& 0.909& 0.579&1.214 & 0.381 & 0.664\\
 950nm & 0.534& 1.041 &0.796& 1.402& 0.628& 0.928
 \end{tabular*}
 {\rule{\temptablewidth}{1pt}}
 \end{center}
 \label{tab:pmc_opti_para}
 \end{table}

\subsubsection{Computation settings}
The computation was distributively performed on a supercomputer BlueGene (CCNI at RPI) with 16,384 nodes (dual CPU at 700MHz). All simulations were running on 1,024 nodes with $10^9$ photons per source, which is the optimal number of photons considering time cost and adequate statistical reliability to compute time-resolved Jacobians in mouse model (Fig. \ref{fig:pMC-recon-stat}). Photon profiles were computed over a 1ns time period. Time gates were simulated with a 200ps width and 20ps gate delay, leading to overall 50 gates.

The inverse calculations were executed on a personal computer (i7 3GHz, 8Gb of RAM), using a Conjugate Gradient algorithm under Matlab. Due to the memory limitation, only three gates were employed for reconstructions. From the TPSFs, the gate with the maximum photon count was selected, as well as the gates having photon count of half the maximum value before and after the maximum gate (Fig. \ref{fig:pmc-jacobian}).

%%-----------------------------------------------------------
\subsubsection{Tomographic reconstruction performance}\label{TG_recon}
As an example, we provide here the time-resolved Jacobian for the center source-detector pair in Fig. \ref{fig:pmc-jacobian}. As expected, the Jacobian shape expands as time increases. Photons collected shortly
after the pulse experience significantly fewer scattering events and probe a small volume, whereas late photons sample larger volumes. This difference in probing volumes provides unique information for reconstruction and is at the origin of the superior performances of
time-gated data types compared to CW data types for DOT as discussed before.
\begin{figure}[htbp]
 \centering
\includegraphics[width=\textwidth]{jacobian}
\caption{Typical Jacobians for three different time gates in a mouse model. First column: TPSF simulated and selected time gates (first half maximum, maximum and last half maximum of TPSFs). Second column: Sagittal plane of the synthetic phantom and normalized Jacobians
corresponding to each selected time gate. Third column: Transverse plane of the murine model with associated normalized Jacobian. The Jacobians are plotted in log scale to provide better visualization.}
\label{fig:pmc-jacobian}
\end{figure}

Functional direct reconstructions were performed on the synthetic mouse model bearing six organs for a variety of scenarios using six wavelengths. Due to the high dynamical range encountered in the different projections computed, we employed a SNR cut-off filter on the measurements selected to compute the inverse problem. Based on heuristic tests, all gates below 500 counts were not considered.

\begin{figure}[htbp]
\centering
\includegraphics[width=\textwidth]{linear_nonlinear}
\caption{Anatomically guided nonlinear sampling: (a) evenly spanned source (black dots) and detectors (red dots) configuration and reconstructed (c) blood volume and (d) oxygen saturation; (b) dense source and detector pairs around the overlapping organs and the reconstructed (e) blood volume and (f) oxygen saturation. The heart and stomach have been delineated with white lines.}
\label{recon_nonlinear}
\end{figure}

An example of functional reconstructions for the case of linear and non-liner source-detector sampling with \emph{a priori} constraints is provided in Fig. \ref{recon_nonlinear}. We investigated the effects of data type (CW or TG), source-detector sampling and soft prior constraints on the accuracy of DOT to retrieve the quantitative bio-distribution of the main functional parameters simulated. The parameters that are reported in this work are the blood volume (BV), the relative oxygen saturation (StO$_2$) and the water content (H$_2$O). The relative errors associated with these parameters are defined by:
\begin{equation}
error_{\overline{BV} }= (\widehat{\overline{BV}}-BV)/BV, \qquad
error_{\overline{H_2O} }= (\widehat{\overline{H_2O}}-H_2O)/H_2O.
\end{equation}
As the StO2 is already a percentage, we provide here the absolute
error:
\begin{equation}
error_{StO_2} = \widehat{\overline{StO_2}}-\overline{StO_2}).
\end{equation}
The standard deviation is calculated as
\begin{equation}
\sigma = \sqrt{\frac{1}{N}\|\hat{X}-\overline{\hat{X}}\|^2},
\end{equation}
where $\hat{X}$ is the reconstructed image for any of the reconstructed quantities. The error and standard deviation functions were computed over the reference volume of each organ derived from the mouse atlas. The values of these error functions for all cases
investigated are provided in Table \ref{tab:result_value}.

The reconstruction using time-gated data with spatial \emph{a priori} presents higher accuracy than using CW data and using time-gated data without spacial constraint. When comparing unconstrained CW and TD reconstruction, one should take into consideration that the introduction of time-domain data is providing better mean value quantification while maintaining the same level of noise. This is in good agreement with results in Section \ref{sec:pmc-single}. The means of \emph{a priori} constrained result are similar to the unconstrained ones while the deviations differs to $~$20\%, indicating that the presence of \emph{a priori} information does improve the noise elimination of the reconstruction.

In all cases investigated herein, the reconstruction using time-gated data with anatomical \emph{a priori} presents overall higher accuracy. Especially, combining TG and spatial constraints provides better estimates than unconstrained CW or TG reconstructions (even though spectral \emph{a priori} was used). Without spatial priors, TG data outperforms CW data, leading to more accurate mean estimates, but suffers from large quantitative spatial deviations within the volume of the organs. This is in good agreement with previous results reported in \cite{Rasmussen2007}.

A significant improvement in the quantitative estimation of the blood volume and the saturation was observed when non-linear dense source patterns were employed. As seen in Fig. \ref{recon_nonlinear}, even source distributions over the whole body lead to under-sampling of important small heterogeneous volume. In this case, the non-linear pattern led to accurate parameters estimation (within 6\% for each organ) whereas linearly distributed patterns did not resolve the heart or the stomach despite the soft spatial constraint.

\begin{sidewaystable}[htbp]
 \tabcolsep 0pt \caption{Errors and standard deviations of the reconstructed maps: a crossed cell in n-Lin or Baye. indicates reconstructions respectively using non-linear source detector sampling or Bayesian formulation in the inverse problem, conversely to linear sampling and unconstrained reconstructions if not marked. All the values herein are provided in percentage (\%).}
 \label{tab:result_value}
\begin{center}
\def\temptablewidth{\textwidth}
{\rule{\temptablewidth}{1pt}}
\begin{tabular*}{\temptablewidth}{@{\extracolsep{\fill}}cccccccccc}
 & Data& Lin.SD &Baye. & Backgd & Heart& liver & kidneys & Stomach
 &lungs\\
\hline
 \multirow{5}{*}{BV}& \multirow{2}{*}{CW}
 &&&17.2$\pm$35.7 & -25.3$\pm$14.4
 & 24.4$\pm$40.5 &-32.0$\pm$20.6 &21.2$\pm$21.5 & 3.1$\pm$27.5\\

 & && $\times$&12.4$\pm$3.4 & -26.2$\pm$0.8& 28.3$\pm$3.1& -36.0$\pm$1.6
 & 20.1$\pm$3.3 & -4.0$\pm$1.4\\

 \cline{2-10}
 &\multirow{3}{*}{TD}
 &&& 15.7$\pm$32.3 & -20.5$\pm$15.2& 22.7$\pm$38.2& -18.2$\pm$19.3
 &12.8$\pm$28.1 & 6.0$\pm$22.4\\

 &&&$\times$& 11.9$\pm$5.9 & -23.8$\pm$1.6&
 25.9$\pm$5.0&-18.6$\pm$2.8&14.6$\pm$4.7&3.9 $\pm$ 3.0\\

 &&$\times$&$\times$&2.3$\pm$5.4 & 3.6$\pm$1.7 & 10.6$\pm$4.0&
 -4.2$\pm$3.0 & 3.0$\pm$4.7 & 13.5$\pm$3.7\\

\hline
 \multirow{5}{*}{StO$_2$}&\multirow{2}{*}{CW}
 &&&3.6$\pm$22.4 & -7.5$\pm$4.0 &28.1$\pm$13.8& 0.9$\pm$70.9&
 4.1$\pm$11.1& -8.4$\pm$6.0\\

 &&&$\times$&-3.0$\pm$0.4 & -8.8$\pm$0.1 &
 14.8$\pm$0.7 &-8.8$\pm$0.5& 11.1$\pm$1.2&
 -5.7$\pm$0.2\\

 \cline{2-10}

 &\multirow{3}{*}{TG}&&&2.9$\pm$16.5& -4.3$\pm$5.3 & 15.2$\pm$10.1& 5.3$\pm$58.3&
 2.8$\pm$11.9 &-10.5$\pm$5.0\\

 &&&$\times$ & 1.8$\pm$0.5 & -4.2$\pm$0.3& 9.80$\pm$1.0
 &-4.3$\pm$0.8 &6.9$\pm$1.0 &-7.6$\pm$0.3 \\

 &&$\times$&$\times$&0.5$\pm$5.4& -1.1$\pm$1.7 &10.7$\pm$4.0
 &-1.8$\pm$3.0 &0.5$\pm$4.7& 1.3$\pm$3.7\\

\hline
 \multirow{5}{*}{H$_2$O}&\multirow{2}{*}{CW}&&&-48.8$\pm$16.7& -58.0 $\pm$23.8& -64.2$\pm$18.2& -48.4$\pm$44.3&-50.1$\pm$21.4& -68.8$\pm$22.8\\

 &&&$\times$&-7.0$\pm$4.8 & -21.8$\pm$ 10.0 &
 -32.8$\pm$7.3 &-27.8$\pm$5.4& -19.1$\pm$12.1&
 -25.7$\pm$2.5\\

 \cline{2-10}

 &\multirow{3}{*}{TG}&&& -49.3$\pm$14.3& -49.1 $\pm$ 9.2 & -60.8$\pm$ 10.2
 & -39.5$\pm$ 10.2& -44.1 $\pm$16.0& -38.9$\pm$ 5.5\\

 &&&$\times$&-17.0$\pm$2.8 & -22.4$\pm$3.2& -30.8$\pm$1.1
 &-35.0$\pm$8.2
 &-10.8$\pm$9.6 &-27.6$\pm$0.3 \\

 &&$\times$&$\times$&-20.9$\pm$3.6& -25.7$\pm$1.3 &-33.8$\pm$6.6
 &-32.4$\pm$5.0 &-5.9$\pm$4.7& -13.7$\pm$6.7
 \end{tabular*}
  {\rule{\temptablewidth}{1pt}}
 \end{center}
\end{sidewaystable}

\subsubsection{Discussion}
In this section, simulations were performed to asses the potential of time-resolved functional DOT for whole-body small animal imaging. To tackle the challenges of small animal imaging, we explored \emph{in silico} the feasibility of combining Monte Carlo based reconstructions, time-gated data and soft \emph{a priori} constraints to quantitatively image the endogenous functional properties of murine model organs. This research presents computational evidence that the combination of time-gated pMC and spatial \emph{a priori} significantly improves the accuracy of diffuse optical tomography in preclinical imaging.

First, the perturbative Monte Carlo formulation is suitable and feasible for functional DOT in preclinical scenarios, where the specimens present a small optically complex model with multiple close organs of significantly low-scattering (bladder
and lungs) or highly absorbing tissues. This method reconstructed wide range of coefficients in tissue accurately and provided satisfactory results synthetically. While only functional parameters related to absorption coefficients were reconstructed in this work, perturbations of structural parameters related to scattering coefficient can be addressed \cite{Li2007, Corlu2003, Corlu2005}.

Second, incorporating \emph{a priori} anatomical constraints provides stable reconstructions and improves quantification. The 3D optical parameters estimation suffers from low resolution, partial volume effect and inter-parameter cross-talk due to the ill-posed nature of DOT. Casting spectral and structural priors into a Bayesian framework yields more accurate and stable calculations to provide enhanced in-vivo functional estimates. Without anatomical guidance, DOT is unable to accurately recover the functional properties of all the organs probed despite spectral priors and dense spatial and temporal datasets.

Third, as a new data type for functional imaging, time gates provide a unique set of information to perform quantitative DOT. Thanks to the different volumes probed by different time gates, time-gated DOT offers superior quantification and resolution compared to CW reconstructions. Time-resolved technology is a promising modality with strength relying on its ability to estimate absolute characterization of tissue optical properties as well as the structural information more precisely. Especially, time-gated technology permits one to acquire both spatially and temporally dense datasets for high-fidelity 3D imaging.

Fourth, this \textit{in silico} investigation provides the first insight on the importance of anatomically guided probe settings. The reconstruction of high absorbing and overlapping tissues is difficult in transmittance settings due to the poor contrast information obtained by non-optimized source-detector configurations. Therefore, image reconstruction from linear source locations suffers not only from fewer photon counts but also from limited sampling of the different organs to be imaged. Even though engineering optimization methods such as the SVA analysis suggest that linear sampling provides a better conditioned inverse problem, non-linear probe setting is of paramount importance to imaging organs such as the heart. However, the punctual probes employed herein suffers from low acquisition time, small penetration depth and low SNR experimentally, thus new strategies for illumination/detection in imaging instrumentation design is in demand. The widefield illumination technique serves as a good candidate to overcome these drawbacks.

\section{Widefield illumination studies}
As stated in Section \ref{sec:intro-chal-wide}, the widefield illumination technique has the potential to be applied to DOT for fast deep tissue imaging. In this section, we develop the MC-based method for the patterned light illumination and detection strategies, and employ this method to investigate the performance of DOT based on time-domain data acquired with simple widefield patterns for the first time. The rationale for this investigation is that the benefits of time-domain tomography observed in point illumination/detection scheme may not readily translate to patterned configuration, as the temporal information may be reduced due to the spatial extent of excitation and detection. Moreover, applying pattern detection strategies to the time domain has not been investigated yet. We show both \emph{in silico} and experimentally that the Monte Carlo-based method as well as the time-gated data with pattern illumination is valid for reconstructions.

\subsection{The Monte Carlo method for modeling widefield illumination}
In the computational aspect, the widefield illumination can be obtained by spatially convolving the punctual forward models with a wider illumination employed. However, this spatial convolution can only be employed when the geometry is rather symmetric, thus it is not applicable in preclinical scenario due to the complex model encountered. The flexible nature of Monte Carlo methods makes it ideal for modeling widefield illumination in complex geometry, and our implementation is elaborated as follow.

A 2D grid having the same resolution as the voxelized geometry is employed to represent the illumination patterns, with each pixel assigned the illumination intensity. Uniformly distributed photons over the illumination area are launched, with the initial weights set to equal the pixel intensity where their initial position falls in. The pattern detection for pMC is achieved by pre-setting the detection patterns in 2D, similar to that for the illumination patterns, and collect the photons exiting the media through the detection patterns.

In the next sections, the pMC method with this widefield illumination/detection capability is employed for tomographic reconstructions. We first establish the merit of the pattern illumination approach in time-domain in simulations and compare the performances of point and patterned light strategies quantitatively. Second, we experimentally demonstrate the potential of widefield structured light illumination and detection scheme for quantitative 3D imaging of absorptive heterogeneities in a murine-like geometry.

\subsection{\emph{in silico} studies}
For the \emph{in silico} studies, a complex absorptive inclusion was simulated within a homogenous medium of size 32mm $\times$32mm and thickness of 20mm (Fig. \ref{fig:ab_wide_setting}a). The background optical properties were set to replicate the experimental settings described in the next paragraph. The absorption heterogeneity had a contrast of 2. The complex absorption structure was positioned in the center of the volume with a thickness of 4mm.

The investigations employed a basic pattern scheme based on a previous study. The patterns represent a large vertical or horizontal bar sliding over the volume of the imaged area (1/2 of the surface is illuminated at any pattern). It may be noted that these patterns do not replicate Fourier modes that are generally sinusoidal in nature. Such simple patterns are easier to implement experimentally, owing to faster switching time due to the absence of gray levels and moreover can be easily applied to complex geometries.

Three experimental strategies were simulated: point source-point detector (po+po), pattern illumination-point detector (pa+po) and pattern illumination-pattern detection (pa+pa). A 2$\times$2cm area (relevant to small animal imaging) was considered for excitation and detection with evenly spanning point sources-detectors pairs. In all \emph{in silico} cases herein, the total number of excitation sources/illumination patterns and point detectors/detection patterns were set at 36 leading to 1,296 combinations (po+po: 36 point source $\times$ 36 point detectors, pa+po: 36 excitation patterns $\times$ 36 point detectors, and pa+po: 36 excitation $\times$ 36 detection patterns).

\begin{figure}[htbp]
\centering
\includegraphics{pmc-4}
\caption{(a) \emph{in silico} phantom with complex absorptive structure. (b) The bar-shape illumination/detection patterns employed in the simulation study sliding along $x$-axis (left column) and $y$-axis (right column).}
\label{fig:ab_wide_setting}
\end{figure}

The computational time for Monte Carlo simulations using the illumination pattern remains generally the same as that using the punctual sources, while a significant increase in computational time
is found for the pattern detections. This increase mainly results from the fact that the probability of a photon being detected is increased due to the wider detection area. Moreover, one photon can be detected by multiple patterns when the patterns are largely overlapped. Both of this will increase the communication time for the slave nodes to send the photons path information to the master nodes.

Examples of typical detector readings for the homogeneous simulations are provided in Fig. \ref{fig:jacobian_wide}. They correspond to the central excitation source/pattern and the central detector/pattern respectively. The Jacobians in the normalized Born scheme associated with these detector readings for the three time gates considered are provided in the same figure. In this study, we selected the maximum, half-rising and half-decaying gates. The gate selection was performed independently for each pattern-detector pair on the homogeneous forward simulation. Thus, for instance, two gates associated with the maximum of the curve can occur at different absolute times for different pattern-detector pairs.

\begin{figure}[htbp]
\centering
\includegraphics[width=\textwidth]{jacobian_wide}
\caption{Simulated detectors readings for the central source-pattern and associated Jacobians for the three gates selected. (a) point source-point detector, (b) patterned illumination-point detector and c) patterned illumination-patterned detection.}
\label{fig:jacobian_wide}
\end{figure}

In all cases presented, the detector readings display the marked temporal features of Temporal Point Spread Functions (TPSF). The use of patterned excitation and/or detection leads to a temporal broadening of the TPSFs (TPSF FWHM: po+po = 0.391ns; pa+po = 0.416ns; pa+pa = 0.415ns) but do not smear the detected photons time of flight statistics significantly. The broadening of the temporal information, associated with the larger volume probed by the detected photons in this configuration, is clearly depicted in the Jacobian spatial distribution. As seen in Fig. \ref{fig:jacobian_wide}, the Jacobians associated with early gates (half-rising maximum) are sensitive to large volumes when widefield structured light is employed, to the extent that few differences can be identified by direct visual comparison of the pa+pa Jacobians for the three gates investigated. Fig. \ref{fig:ab_wide_recon} shows a comparison of the optical reconstructions based on these time resolved datasets. For a similar dataset size, pa+po strategy outperforms po+po both in quantification and resolution. Whereas po+po and pa+po offer relatively similar performances. Even though, pa+po is characterized by spatially extended Jacobians, the optical reconstructions based on this strategy allows for quantitative three dimensional imaging as depicted in Fig. \ref{fig:ab_wide_recon}c.

\begin{figure}[htbp]
\centering
\includegraphics{ab_wide_recon}
\caption{Optical reconstructions of the phantom for time-gated data with (a) point source-point detector, (b) patterned excitation-point detector and (c) patterned excitation - patterned detector settings. The isovolume was set at 40\% of the maximum reconstructed value. The projections of the central reconstructed plane are provided on each boundary for each case.}
\label{fig:ab_wide_recon}
\end{figure}

\subsection{Experimental Studies}
\subsubsection{\label{sec:wide-instrument}System design for time-resolved widefield tomography}
The algorithm was experimentally validated using a time-domain small animal molecular imaging system, which is schematically depicted in Fig.~\ref{fig:ab_exp}. The system employs a tunable Ti-Sapphire laser (Mai Tai HP, Spectra-Physics, CA, USA) which generates 100fs pulses at 80MHz in the NIR spectral band (690nm-1,020nm) as the source. The average laser power has dynamic range of 760mW-3W over its tuning range. A computer-controlled power controller (Application Note 30, Newport, CA, USA) was used in conjunction with the laser to maintain the incident laser power at a constant value with an accuracy of 100mW when tuning the laser during the imaging session. The source was injected into a 400mm multi-mode fiber with 0.22 numerical aperture which guides the pulses to a 97/3 non-polarizing beam-splitter (Ozoptics, ON, Canada). The 3\% channel was connected to a variable attenuator (Ocean Optics, FL, USA) and was focused directly on a diffusing paper on the imaging stage through a 400$\mu$m multi-mode fiber. This channel was used as a temporal reference measured simultaneously during acquisition. The variable attenuator was used to adjust the signal level on the reference channel while tuning to correct for the change in spectral sensitivity of the detector. This allows for an in-experiment calibration of the incident power and $t_0$ of the impulse response function (IRF), reducing the estimation errors due to uncertainties in the time of injection [35]. The 97\% channel was injected into 15x beam expander (BE15M-B, Thorlabs). The expanded beam was incident on digital micro-mirror device based light processing (DLP) board (Discovery 1100, Texas Instruments). The DLP comprises of 1024$\times$768 micro-mirrors which were independently controlled by pulse-width modulation to generate a maximum of 256 grayscale levels. The pixelized nature of the micro-mirrors however resulted in 2D diffraction patterns in the reflected beam where the higher order diffraction beams appear as ghost images [36]. In this system, the higher-orders were re-imaged on the imaging chamber using a bi-convex lens (75mm focal length) to provide an area of excitation spanning 40mm$\times$25mm. It should be noted that the position of the lens used can be modified to optimize the area of illumination on the imaging chamber. Moreover, as the tip of the fiber is directly imaged on the chamber, the intensity profile of the excitation pattern was non-uniform which necessitates the calibration of the excitation pattern within each experiment protocol. Furthermore, the pattern obtained can thus be supplied directly to the Monte Carlo model allowing an accurate modeling of the experimental settings.

The time-gated detection system was implemented using an ultrafast gated intensified CCD camera (Picostar HR, LaVision GmbH, Germany) as the detector. The synchronization of the laser pulses with the intensifier shutter was achieved using an optical trigger (OCF-401, Becker \& Hickl GmbH, Germany). The pulse train was conditioned using the HRI-delay unit (Kentech Instruments, UK) to potentially achieve 1ps resolution in the temporal measurements over a 50ns scan range. The output from the delay unit was connected to the High-Rate Imager (HRI) which controled the intensifier gating and gain modulation. Intensifier gating refers to the operation of the shutter where the photons are integrated at each time bin over a specified duration, referred to as the gate-width. The gate-width can be set to a preset value ranging from 200ps to 1,000ps, (COMB modes) or controlled through external signals (RF or logical triggering). In this implementation the HRI was operated in the COMB mode where the gate was shifted over the measurement window at specific intervals to measure the time-gated data. The measured signals were further amplified by changing the gain voltage across the microchannel plate (MCP) in the intensifier before impinging on a P43 phosphor screen. The images formed on the screen were imaged by the 12-bit CCD at a resolution of 1,376 $\times$ 1,040. The ICCD camera can detect a maximum of 4,096 photons at each gate and the recorded images were binned over 8 $\times$ 8 pixels post-acquisition to create 1mm$^2$ detector measurements.

\begin{figure}[htbp]
\centering
\includegraphics{pmc-1}
\caption{Schematic of the optical tomography system.}
\label{fig:ab_exp}
\end{figure}

\subsubsection{Reconstructing absorptive perturbation}
A polycarbonate tank of dimensions 80mm$\times$50mm$\times$20mm was filled with a mixture of Intralipid-20\% (Sigma-Aldrich, USA) diluted in water and India Ink to simulate bulk murine tissue optical parameters ($\mu_a$=0.005mm$^{-1}$ and $\mu_s'$=1.16mm$^{-1}$). The phantom hold 3 cylindrical glass inclusions having a maximum diameter of 7mm distributed over the volume as depicted in Fig. \ref{fig:wide_phantom}b. The inclusions were filled with the same mixture as the background to obtain homogeneous measurements and a liquid absorber was added to create a contrast of 8x for objects II and III, and 4x for object I for heterogeneous measurements. The phantom was illuminated over a 45mm by 35mm surface with the widefield source created using the platform described above at 700nm. The transmitted modulated time-resolved light was detected with an ultra-fast gated CCD camera. The acquisition parameters of the camera were set to 300ps gates at 40ps intervals with a high voltage of 510V. The overall time of acquisition for 36 illumination patterns was less than 12min.

\begin{figure}[htbp]
\centering
\includegraphics{pmc-5}
\caption{(a) Experimental phantom with three absorptive inclusions. The red marking corresponds to the illuminated surface and point detector distribution. The white dot lines delineate the reconstructed volume. (b) Sample illumination patterns in the experiment sliding along $x$-axis (left column) and $y$-axis (right column).}
\label{fig:wide_phantom}
\end{figure}

The 3D experimental reconstructions for absorption heterogeneities is shown in Fig. \ref{fig:ab_wide_exp_recon}. In this investigation we were limited to widefield strategies as it was not possible to concurrently acquire point excitation and widefield excitation data with our current instrument. Note also that the patterned detection scheme was obtained from the same dataset by applying integrating masks on the CCD camera output. However, widefield pattern detection can be performed with one single detector associated with a second digital micromirror device (DMD) in the detection optical chain \cite{Belanger2010}. The benefits are then on the instrumental side: cost, use of more sensitive detector (PMT), potential for spectral multiplexing, etc..

\begin{figure}[htbp]
\centering
\includegraphics{ab_wide_exp_recon}
\caption{Optical reconstructions of the phantom for the pattern source-point detector strategy. (a) CW data, (b) time-gated data. The isovolume was set at 50\% of the maximum reconstructed value. The ray-sum projections are provided on the sides of the reconstructed volume.}
\label{fig:ab_wide_exp_recon}
\end{figure}

Table \ref{tab:ab_wide_exp} summarizes a quantitative assessment of the different experimental strategies investigated. We computed two metrics based on the 3D reconstruction. First, a quantitative accuracy metric calculated as the error in mean absorption coefficient reconstructed over the 50\% isovolume versus the known value was computed. Second, a resolution metric defined as the estimated radius of the cylindrical object reconstructed compared to the true radius of the object (7mm) was estimated.

\begin{table}[htbp]
 \tabcolsep 0pt \caption{\label{tab:ab_wide_exp}Quantitative and resolution errors in the 3D reconstructions compared to the known values of the experimental heterogeneities. pa: pattern; po: point.} \vspace*{-12pt}
\begin{center}
\def\temptablewidth{\textwidth}
{\rule{\temptablewidth}{1pt}}
\begin{tabular*}{\temptablewidth}{@{\extracolsep{\fill}}ccccccccc}
 \multirow{2}{*}{Source}&\multirow{2}{*}{Detector} &\multirow{2}{*}{Datatype} &\multicolumn{3}{c}{Estimation Error}& \multicolumn{3}{c}{Diameter (mm)}\\
 & &&Incl I &Incl II &Incl III &Incl I &Incl II & Incl III \\
\hline
 \multirow{2}{*}{pa}&\multirow{2}{*}{po} &TG &-30.3 & 3.6 & 2.7 &6.6 &7.8 & 9.1 \\
 & &CW &17.7 & 5.3 & 26.1 & 8.0 & 8.4 &14.1 \\
 \multirow{2}{*}{pa}&\multirow{2}{*}{pa} &TG &6.5 & 42.4 & 17.5 &17.1 &14.8 & 14.1 \\
 & &CW &31.9 & 57.7 & 32.9 & 14.6 & 16.6 &11.4
 \end{tabular*}
{\rule{\temptablewidth}{1pt}}
 \end{center}
 \end{table}

In all the cases investigated, the three objects were resolved in reconstruction. The pa+po with TG datasets provided more accurate reconstructions than their CW equivalents. Thus even with extended illumination/detection strategies, the additional information provided by TG data types was retained and allowed for higher fidelity reconstructions than CW data types. In particular, the improvement in resolving the depth information is clearly demonstrated in the reconstruction of object III.

As expected from \emph{in silico} simulations, the experimental reconstructions based on illumination and detection patterns were characterized by a poorer resolution and quantification. Only in the case of object I, pa+po was less quantitatively accurate as object I was over-reconstructed, leading to a smaller reconstructed diameter and hence, a higher mean absorption than expected.

% needs to be rewritten
\subsubsection{Reconstructing absorption and scattering perturbation simultaneously}
In this experiment we investigated the performance of the pMC-based model for tomographic reconstruction of absorption and scattering coefficients using the time-gate data on a murine model phantom. The phantom (Fig. \ref{fig:absc_exp_setting}a and \ref{fig:absc_exp_setting}b) consisted of a 2cm thick polycarbonate tank carrying a mixture of Intralipid-20\% and India ink ($\mu_a$= 0.05cm$^{-1}$; $\mu_s$' = 9cm$^{-1}$). Two tubes with 10mm inner diameter were suspended in the tank with India ink and Intralipid-20\% solution added to tubes 1 and 2 to simulate localized contrast in absorption ($\mu_a$ = 0.4cm$^{-1}$; $\mu_s$' = 9cm$^{-1}$) and scattering ($\mu_a$ = 0.05cm$^{-1}$; $\mu_s$' = 18cm$^{-1}$) respectively. The inclusions were selected to simulate perturbations comparable to the size of large organs (e.g. liver, lungs) in murine models.

\begin{figure}[htbp]
\centering
\includegraphics[width=\textwidth]{pmc-2}
\caption{Reconstruction of $\mu_a$ and $\mu_s$ \emph{in vitro}: (a) Phantom design. (b) Area of pattern illumination. (c)-(d) TPSF measured at the green and cyan detectors respectively. (e)-(g) Normalized Born contrast measures at TG1, TG2 and TG3 respectively.}
\label{fig:absc_exp_setting}
\end{figure}

A set of basic bar shaped patterns spanning half of the excited surface (40mm $\times$ 25mm) were employed as the widefield illumination scheme for tomographic reconstruction (Fig. \ref{fig:absc_exp_setting}a). The experiment protocol employed 36 'bar' patterns (18 patterns along the $x$-axis and 18 patterns along the $y$-axis) as the source. The intensifier was operated with a gatewidth of 300ps at 20ps intervals spanning a 2.2ns time window and two sets of measurements were acquired; first, with the tube containing the background mixture (homogeneous) and second, with the added India ink and Intralipid solution (heterogeneous). The entire acquisition protocol was completed in 9 minutes. The acquired images were postprocessed to generate 1mm $\times$ 1mm detectors. 88 point detectors with a separation of 3mm along the $x$-axis and 5mm along the $y$-axis were selected for the reconstruction.

Fig. \ref{fig:absc_exp_setting}b shows the comparison of the homogenous and heterogeneous TPSFs at two point detectors located over the absorption and scattering perturbation. It should be noted that the absorptive perturbation produces higher contrast for later gates while the scattering perturbation has predominant contrast in the early gates with minimal contrast observed in later gates. The normalized born contrast function at three time gates across the excited volume (TG1: 20\% of peak intensity on the rising edge, TG2: peak intensity and TG3: 50\% of peak intensity on the falling edge) shows a distinct separation of contrast due to the two objects. The scattering object presented predominant contrast during the early gate, while showing minimal contrast during the late gates (Fig. \ref{fig:absc_exp_setting}c).
\begin{figure}[htbp]
\centering
\includegraphics{pmc-3}
\caption{50\% iso-volumes of reconstructed perturbations in optical
properties.}
\label{fig:absc_exp_result}
\end{figure}

To solve the inverse problem, the measurement vector was constructed by heuristically selecting 6 time gates for each detector (8\%, 12\% and 17\% of the peak on the rising edge; peak; 80\% and 60\% of the peak on the falling edge). It should be noted that the time gates were selected based on the homogenous measurement. The phantom volume was modeled as a slab discretized into 2mm $\times$ 2mm $\times$ 2mm voxels and the Jacobians for perturbations in the scattering and absorption coefficients were computed using the average optical properties. The MC simulations using $2 \times 10^{10}$ photons were launched on a supercomputer BlueGene (CCNI at RPI). The Jacobians were scaled by their mean-column value to reduce the inter-parameter cross talk \cite{Intes2005a}. A spatially varying regularization was applied while solving the inverse problem\cite{Pogue1999}. The regularization parameter was empirically determined to have a value of $10^3$ on the source detector layers and $5\times10^2$ on the penultimate layers. The problem was solved using a least-squares algorithm (lsqr, MATLAB). The 3D visualization of the reconstructed optical properties is shown in Fig.~\ref{fig:absc_exp_result}. The mean value of the 50\% isovolume for both parameters was found to have an error less than 5\% (absorption coefficient, inclusion 1: expected value = 0.35cm$^{-1}$, reconstructed mean value = 0.33 cm$^{-1}$; scattering coefficient, inclusion 2: expected value = 9.0cm$^{-1}$, reconstructed mean value = 8.92cm$^{-1}$).

Besides the DOT performance consideration, widefield strategies potentially offer numerous experimental advantages. First, the spatial distribution of collected photons in transmittance is characterized by a smaller dynamic range compared to a point source excitation scheme. This enables data collection with satisfactory signal-to-noise ratio for stable reconstruction over a extended spatial area compare to point illumination strategies. Second, the strategies employing patterned light are less sensitive to the boundaries than discrete optodes as demonstrated in the Jacobian spatial distribution of Fig. \ref{fig:jacobian_wide}. Thus, structured patterned strategies are potentially more appropriate for monitoring deep tissue without interference from superficial layers, benefitting applications such as brain functional monitoring. In addition, the time of acquisition for spatial information over large volumes is considerably reduced compared to raster scan or fiber coupled techniques. So far the technique is based on time-gated data and for absorption heterogeneities, but the technique can be extended to any time-resolved datasets and contrast functions.

Lastly, the modeling of illumination patterns is achieved thanks to the Monte Carlo method's excellent flexibility. In this section, the simple half-space bar shaped patterns were employed, however, arbitrary patterns can be simulated using this approach, which opens the door for pattern optimization on complex models. Moreover, the pa+po setting provides better performance of reconstruction in both resolution and quantification than the pa+pa setting, while using the pa+po setting will not result in extended computational time for the Monte Carlo simulations. Thus the computational efficiency of this method is retained while giving significant experimental benefits.

\section{Summary}
In this chapter, a Monte-Carlo-based time-resolved approach was developed for accurate reconstruction of functional parameters. The low time efficiency limitation of Monte Carlo simulations was resolved by the use of massive parallel computing techniques. This method was applied to simulated data for imaging optical properties in a synthetic mouse model. Then we explored the feasibility of employing soft \textit{a priori} constraints to guide (a) the probe configuration, and (b) the reconstruction iteration, to quantitatively image the endogenous functional properties of \textit{in silico} murine model organs. The combination of time-gated pMC and spatial \emph{a priori} significantly improved the accuracy of DOT in preclinical scenarios where the diffusion theory fails. The proposed method was then extended to widefield illumination settings for experimental validation. Optical properties were successfully reconstructed using the data acquired by an instrumentation platform with widefield illumination. This reconstruction strategy can be used to produce an \emph{in vivo} functional atlas of murine model, which will be extremely valuable for the biomedical community. Moreover, the proposed method can be extended to fluorescence molecular imaging, which is discussed in the next chapter. 